library(foreign); library(survival); library(ggplot2); library(coefplot); library(ggpubr)
library(mediation); library(MASS); library(survey); library(margins); library(effects); library(colorRamps)
library(readstata13); library(reshape2); library(lattice); library(scales); library(ggthemes); library(RColorBrewer)
library(effects);library(readstata13); library(stargazer); library(gridExtra)
####read in data
setwd('~<insert working directory here>')
data<- read.dta13(file='shepperd-vonstein-final.dta', convert.underscore=TRUE)

###generate variable indicating treatment group###
data$treatmentgroupfactor <- factor(data$treatmentgroup, levels = c("control", "il", "moral", "reputational"), 
  labels = c("Control", "Intl. Law", "Moral", "Reputational"))
##############generate Figure 1a#################
data$female1 <-factor(data$female, levels=c("0", "1"),labels =c("Male", "Female"))
female_bar <- ggplot(data, aes(x = treatmentgroup, group=female1)) +
  geom_bar(aes(y = ..prop.., fill = factor(..x..)), stat="count") +
  facet_grid(~ female1) +
  scale_fill_brewer(palette = "Set1", labels=c("Control", "Intl. law", "Moral", "Reputational")) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), 
        legend.position = "bottom", axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank(), axis.line = element_blank()) + 
  labs(y="Proportion in Group") + guides(fill=guide_legend(title="Treatment Group"))
female_bar

data$agegroup1 <-factor(data$agegroup, levels=c("1", "2", "3", "4"),labels =c("18-34", "35-49", "50-64", "65+"))
age_bar <- ggplot(data, aes(x = treatmentgroup, group=agegroup1)) +
  geom_bar(aes(y = ..prop.., fill = factor(..x..)), stat="count") +
  facet_grid(~ agegroup1) +
  scale_fill_brewer(palette = "Set1", labels=c("Control", "Intl. law", "Moral", "Reputational")) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), 
        legend.position = "bottom", axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank(), axis.line = element_blank()) + 
  labs(y="Proportion in Group") + guides(fill=guide_legend(title="Treatment Group"))
age_bar

###combine both plots to make Figure 1a###########
ggarrange(female_bar, age_bar, ncol = 2, nrow = 1,widths = c(1, 2.5) )

##############generate Figure 2a#################
data<- read.dta13(file='shepperd-vonstein-final.dta', convert.underscore=TRUE)
#subsetting to leave out those who responded "prefer not to say" (n = 44)
data<- subset(data, education1>=1&education1<=3)
data$educationfactor <-factor(data$education1, levels=c("1", "2", "3"),
 labels =c("Unskilled Trade or Less", "Skilled Trade School", "University"))
education_bar <- ggplot(data, aes(x = treatmentgroup, group=educationfactor)) +
  geom_bar(aes(y = ..prop.., fill = factor(..x..)), stat="count") +
  facet_grid(~ educationfactor) +
  scale_fill_brewer(palette = "Set1", labels=c("Control", "Intl. law", "Moral", "Reputational")) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), 
  axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank(), 
  axis.line = element_blank()) + theme(legend.position="none")+
  labs(y="Proportion in Group") + guides(fill=guide_legend(title="Treatment Group"))
education_bar

####################generate Figure 3a#############
data<- read.dta13(file='shepperd-vonstein-final.dta', convert.underscore=TRUE)
data$incomefactor <- factor(data$income, levels = c("1", "2", "3", "4","5"),
  labels = c("<$24,999", "$25,000-$49,999", "$50,000-$74,999", "$75,000-$99,999", "$100,00+"))
income_bar <- ggplot(data, aes(x = treatmentgroup, group=incomefactor)) +
  geom_bar(aes(y = ..prop.., fill = factor(..x..)), stat="count") +
  facet_grid(~ incomefactor) +
  scale_fill_brewer(palette = "Set1", labels=c("Control", "Intl. law", "Moral", "Reputational")) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), legend.position = "top", axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank(), axis.line = element_blank()) + 
  labs(y="") + theme(legend.position="none")
guides(fill=guide_legend(title="Treatment Group"))
income_bar
####################generate Figure 4a#############
data<- read.dta13(file='shepperd-vonstein-final.dta', convert.underscore=TRUE)
data$ethnicdist <- factor(data$ethnicdistcollapsed, levels = c("0", "1", "2"),
 labels = c("Strongly/disagree", "Neutral", "Strongly/agree"))
ethnic_bar <- ggplot(data, aes(x = treatmentgroup, group=ethnicdist)) +
  geom_bar(aes(y = ..prop.., fill = factor(..x..)), stat="count") + facet_grid(~ ethnicdist) + 
  scale_fill_brewer(palette = "Set1", labels=c("Control", "Intl. law", "Moral", "Reputational")) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), 
  legend.position = "none", axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank(), 
  axis.line = element_blank()) + labs(y="Proportion in Group") 
ethnic_bar

####################generate Figure 5a#############
data<- read.dta13(file='shepperd-vonstein-final.dta', convert.underscore=TRUE)
political_bar <- ggplot(data, aes(x = treatmentgroup, group=political)) +
  geom_bar(aes(y = ..prop.., fill = factor(..x..)), stat="count") +
  facet_grid(~ political) + scale_fill_brewer(palette = "Set1", name = "Treatment Group", labels=c("Control", "Intl. law", 
  "Moral", "Reputational")) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), 
        legend.position = "none", axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.x=element_blank(), axis.line = element_blank()) + 
  labs(y="Proportion in Group")
political_bar

####################generate Table 10a#############
petition.base<-glm(petition1~ilframe+moralframe+reputationalframe, data=data, family=binomial(link="probit"));summary(petition.base)
petition.extended<-glm(petition1~ilframe+moralframe+reputationalframe+education1+income+agegroup+ethnicdist+female+labor+greens+
    otherindependent+ dontknowwontvote, data=data, family=binomial(link="probit"));summary(petition.extended)
 
protest.base<-glm(protest1~ilframe+moralframe+reputationalframe, data=data, family=binomial(link="probit"));summary(protest.base)
protest.extended<-glm(protest1~ilframe+moralframe+reputationalframe+education1+income+agegroup+ethnicdist+female+labor+greens+
   otherindependent+ dontknowwontvote, data=data, family=binomial(link="probit"));summary(protest.extended)

donate.base<-glm(donate1~ilframe+moralframe+reputationalframe, data=data, family=binomial(link="probit"));summary(donate.base)
donate.extended<-glm(donate1~ilframe+moralframe+reputationalframe+education1+income+agegroup+ethnicdist+female+labor+greens+
  otherindependent+ dontknowwontvote, data=data, family=binomial(link="probit"));summary(donate.extended)

stargazer(petition.base, petition.extended, protest.base, protest.extended, donate.base, donate.extended, ci=TRUE, type="html",out="action2.doc",star.cutoffs = c(.05, .01))

###Table 11a, Wald tests
###for petition###
linearHypothesis(petition.base, "ilframe=moralframe")
linearHypothesis(petition.base, "ilframe=reputationalframe")
linearHypothesis(petition.base, "moralframe=reputationalframe")
linearHypothesis(petition.extended, "ilframe=moralframe")
linearHypothesis(petition.extended, "ilframe=reputationalframe")
linearHypothesis(petition.extended, "moralframe=reputationalframe")
###for protest###
linearHypothesis(protest.base, "ilframe=moralframe")
linearHypothesis(protest.base, "ilframe=reputationalframe")
linearHypothesis(protest.base, "moralframe=reputationalframe")
linearHypothesis(protest.extended, "ilframe=moralframe")
linearHypothesis(protest.extended, "ilframe=reputationalframe")
linearHypothesis(protest.extended, "moralframe=reputationalframe")
###for donation###
linearHypothesis(donate.base, "ilframe=moralframe")
linearHypothesis(donate.base, "ilframe=reputationalframe")
linearHypothesis(donate.base, "moralframe=reputationalframe")
linearHypothesis(donate.extended, "ilframe=moralframe")
linearHypothesis(donate.extended, "ilframe=reputationalframe")
linearHypothesis(donate.extended, "moralframe=reputationalframe")

###Table 12a, power analysis###
###create subsets###
control.subset<-subset(data,noframe==1)
il.subset<-subset(data,ilframe==1)
moral.subset<-subset(data,moralframe==1)
reputational.subset<-subset(data,reputationalframe==1)
#calculate N for each group
mean.controlpolicyapproval1<-mean(control.subset$policyapproval1)
n.control<-nrow(data[data$noframe == "1",])
n.il<-nrow(data[data$ilframe == "1",])
n.moral<-nrow(data[data$moralframe == "1",])
n.reputational<-nrow(data[data$reputationalframe == "1",])
###calculate mean for each group###
mean.controlpetition<-mean(control.subset$petition1)
mean.ilpetition<-mean(il.subset$petition1)
mean.moralpetition<-mean(moral.subset$petition1)
mean.reputationalpetition<-mean(reputational.subset$petition1)
###Power analysis reported in Table 12a -- petition
pwr.t2n.test(n1 = n.il, n2=n.control , d = (mean.ilpetition-mean.controlpetition) , sig.level =.05 ); summary(n.il+n.control) #IL-control
pwr.t2n.test(n1 = n.il, n2=n.moral , d = (mean.ilpetition-mean.moralpetition) , sig.level =.05 );summary(n.il+n.moral)#IL-moral
pwr.t2n.test(n1 = n.il, n2=n.reputational , d = (mean.ilpetition-mean.reputationalpetition) , sig.level =.05 );summary(n.il+n.reputational)#IL-reputational
pwr.t2n.test(n1 = n.moral, n2=n.control , d = (mean.moralpetition-mean.controlpetition) , sig.level =.05 );summary(n.moral+n.control)#Moral-reputational
pwr.t2n.test(n1 = n.moral, n2=n.reputational , d = (mean.moralpetition-mean.reputationalpetition) , sig.level =.05 );summary(n.moral+n.reputational)#Moral-reputational
pwr.t2n.test(n1 = n.reputational, n2=n.control , d = (mean.reputationalpetition-mean.controlpetition) , sig.level =.05 );summary(n.reputational+n.control)#Reputational-control
###calculate mean for each group###
mean.controlprotest<-mean(control.subset$protest1)
mean.ilprotest<-mean(il.subset$protest1)
mean.moralprotest<-mean(moral.subset$protest1)
mean.reputationalprotest<-mean(reputational.subset$protest1)
###Power analysis reported in Table 12a -- protest
pwr.t2n.test(n1 = n.il, n2=n.control , d = (mean.ilprotest-mean.controlprotest) , sig.level =.05 ); summary(n.il+n.control) #IL-control
pwr.t2n.test(n1 = n.il, n2=n.moral , d = (mean.ilprotest-mean.moralprotest) , sig.level =.05 );summary(n.il+n.moral)#IL-moral
pwr.t2n.test(n1 = n.il, n2=n.reputational , d = (mean.ilprotest-mean.reputationalprotest) , sig.level =.05 );summary(n.il+n.reputational)#IL-reputational
pwr.t2n.test(n1 = n.moral, n2=n.control , d = (mean.moralprotest-mean.controlprotest) , sig.level =.05 );summary(n.moral+n.control)#Moral-reputational
pwr.t2n.test(n1 = n.moral, n2=n.reputational , d = (mean.moralprotest-mean.reputationalprotest) , sig.level =.05 );summary(n.moral+n.reputational)#Moral-reputational
pwr.t2n.test(n1 = n.reputational, n2=n.control , d = (mean.reputationalprotest-mean.controlprotest) , sig.level =.05 );summary(n.reputational+n.control)#Reputational-control
###calculate mean for each group###
mean.controldonate<-mean(control.subset$donate1)
mean.ildonate<-mean(il.subset$donate1)
mean.moraldonate<-mean(moral.subset$donate1)
mean.reputationaldonate<-mean(reputational.subset$donate1)
###Power analysis reported in Table 12a -- donate
pwr.t2n.test(n1 = n.il, n2=n.control , d = (mean.ildonate-mean.controldonate) , sig.level =.05 ); summary(n.il+n.control) #IL-control
pwr.t2n.test(n1 = n.il, n2=n.moral , d = (mean.ildonate-mean.moraldonate) , sig.level =.05 );summary(n.il+n.moral)#IL-moral
pwr.t2n.test(n1 = n.il, n2=n.reputational , d = (mean.ildonate-mean.reputationaldonate) , sig.level =.05 );summary(n.il+n.reputational)#IL-reputational
pwr.t2n.test(n1 = n.moral, n2=n.control , d = (mean.moraldonate-mean.controldonate) , sig.level =.05 );summary(n.moral+n.control)#Moral-reputational
pwr.t2n.test(n1 = n.moral, n2=n.reputational , d = (mean.moraldonate-mean.reputationaldonate) , sig.level =.05 );summary(n.moral+n.reputational)#Moral-reputational
pwr.t2n.test(n1 = n.reputational, n2=n.control , d = (mean.reputationaldonate-mean.controldonate) , sig.level =.05 );summary(n.reputational+n.control)#Reputational-control

####################Figure 13a############
ilcontrol.subset<-subset(data,ilframe==1|noframe==1)
model.m.ilcontrol<- lm(policyapproval1~ilframe, data = ilcontrol.subset)
model.y.ilcontrol<-glm(petition1~ilframe+policyapproval1, data = ilcontrol.subset, family=binomial(link="probit"))
mediatedilcontrol.petition<- mediate(model.m.ilcontrol, model.y.ilcontrol, treat = "ilframe", mediator = "policyapproval1", boot=TRUE, sims=2000,conf.level = .95)
summary(mediatedilcontrol.petition);plot(mediatedilcontrol.petition)

#il group vs. moral
ilmoral.subset<-subset(data,ilframe==1|moralframe==1)
model.m.ilmoral<- lm(policyapproval1~ilframe, data = ilmoral.subset)
model.y.ilmoral<-glm(petition1~ilframe+policyapproval1,  data = ilmoral.subset, family=binomial(link="probit"))
mediatedilmoral.petition<- mediate(model.m.ilmoral, model.y.ilmoral, treat = "ilframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedilmoral.petition); plot(mediatedilmoral.petition)

#il group vs. reputational
ilreputational.subset<-subset(data,ilframe==1|reputationalframe==1)
model.m.ilreputational<- lm(policyapproval1~ilframe, data = ilreputational.subset)
model.y.ilreputational<-glm(petition1~ilframe+policyapproval1, data = ilreputational.subset,family=binomial(link="probit"))
mediatedilreputational.petition<- mediate(model.m.ilreputational, model.y.ilreputational, treat = "ilframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedilreputational.petition); plot(mediatedilreputational.petition)

#moral group vs. control
moralcontrol.subset<-subset(data,moralframe==1|noframe==1)
model.m.moralcontrol<- lm(policyapproval1~moralframe, data = moralcontrol.subset)
model.y.moralcontrol<-glm(petition1~moralframe+policyapproval1,  data = moralcontrol.subset,family=binomial(link="probit"))
mediatedmoralcontrol.petition<- mediate(model.m.moralcontrol, model.y.moralcontrol, treat = "moralframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedmoralcontrol.petition);plot(mediatedmoralcontrol.petition)

#moral group vs. il group: no estimation -- this is just the inverse of il group vs. moral (above)

#moral group vs. reputational
moralreputational.subset<-subset(data,moralframe==1|reputationalframe==1)
model.m.moralreputational<- lm(policyapproval1~moralframe, data = moralreputational.subset)
model.y.moralreputational<-glm(petition1~moralframe+policyapproval1,  data = moralreputational.subset,family=binomial(link="probit"))
mediatedmoralreputational.petition<- mediate(model.m.moralreputational, model.y.moralreputational, treat = "moralframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedmoralreputational.petition); summary(mediatedmoralreputational.petition)

#reputational group vs. control
reputationalcontrol.subset<-subset(data,reputationalframe==1|noframe==1)
model.m.reputationalcontrol<- lm(policyapproval1~reputationalframe, data = reputationalcontrol.subset)
model.y.reputationalcontrol<-glm(petition1~reputationalframe+policyapproval1,  data = reputationalcontrol.subset,family=binomial(link="probit"))
mediatedreputationalcontrol.petition<- mediate(model.m.reputationalcontrol, model.y.reputationalcontrol, treat = "reputationalframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedreputationalcontrol.petition);plot(mediatedreputationalcontrol.petition)

#reputational group vs il group: no estimation -- this is just the inverse of il group vs. reputational (above)
#reputational group vs. moral group: no estimation -- this is just the inverse of moral group vs. reputational (above)

#########Figure 14a############
sigma.x.ilcontrol<-sd(ilcontrol.subset$ilframe)#calculate std dev of x for power analysis
sigma.m.ilcontrol<-sd(ilcontrol.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.ilcontrol<-sd(ilcontrol.subset$petition1)#calculate std dev of y for power analysis
summary(n.il+n.control)#sample size to report
a<- summary(model.m.ilcontrol)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.ilcontrol)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.il+n.control), power = NULL, a =a , b =b, varx =var(ilcontrol.subset$ilframe), vary=var(ilcontrol.subset$petition1), 
             varm =var(ilcontrol.subset$policyapproval1), alpha = .05, interval = NULL)

#il-moral
sigma.x.ilmoral<-sd(ilmoral.subset$ilframe)#calculate std dev of x for power analysis
sigma.m.ilmoral<-sd(ilmoral.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.ilmoral<-sd(ilmoral.subset$petition1)#calculate std dev of y for power analyais
summary(n.il+n.moral)#sample size to report
a<- summary(model.m.ilmoral)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.ilmoral)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.il+n.moral), power = NULL, a =a , b =b, varx =var(ilmoral.subset$ilframe), vary=var(ilmoral.subset$petition1), 
             varm =var(ilmoral.subset$policyapproval1), alpha = .05, interval = NULL)

#il-reputational
sigma.x.ilreputational<-sd(ilreputational.subset$ilframe)#calculate std dev of x for power analysis
sigma.m.ilreputational<-sd(ilreputational.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.ilreputational<-sd(ilreputational.subset$petition1)#calculate std dev of y for power analyais
summary(n.il+n.reputational)#sample size to report
a<- summary(model.m.ilreputational)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.ilreputational)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.il+n.reputational), power = NULL, a =a , b =b, varx =var(ilreputational.subset$ilframe), vary=var(ilreputational.subset$petition1), 
             varm =var(ilreputational.subset$policyapproval1), alpha = .05, interval = NULL)

#moral-control
sigma.x.moralcontrol<-sd(moralcontrol.subset$moralframe)#calculate std dev of x for power analysis
sigma.m.moralcontrol<-sd(moralcontrol.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.moralcontrol<-sd(moralcontrol.subset$petition1)#calculate std dev of y for power analyais
summary(n.moral+n.control)#sample size to report
a<- summary(model.m.moralcontrol)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.moralcontrol)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.moral+n.control), power = NULL, a =a , b =b, varx =var(moralcontrol.subset$moralframe), vary=var(moralcontrol.subset$petition1), 
             varm =var(moralcontrol.subset$policyapproval1), alpha = .05, interval = NULL)

###moral-reputational
sigma.x.moralreputational<-sd(moralreputational.subset$moralframe)#calculate std dev of x for power analysis
sigma.m.moralreputational<-sd(moralreputational.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.moralreputational<-sd(moralreputational.subset$petition1)#calculate std dev of y for power analyais
summary(n.moral+n.reputational)#sample size to report
a<- summary(model.m.moralreputational)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.moralreputational)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.moral+n.reputational), power = NULL, a =a , b =b, varx =var(moralreputational.subset$moralframe), vary=var(moralreputational.subset$petition1), 
             varm =var(moralreputational.subset$policyapproval1), alpha = .05, interval = NULL)

##reputational-control
sigma.x.reputationalcontrol<-sd(reputationalcontrol.subset$reputationalframe)#calculate std dev of x for power analysis
sigma.m.reputationalcontrol<-sd(reputationalcontrol.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.reputationalcontrol<-sd(reputationalcontrol.subset$petition1)#calculate std dev of y for power analyais
summary(n.reputational+n.control)#sample size to report
a<- summary(model.m.reputationalcontrol)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.reputationalcontrol)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.reputational+n.control), power = NULL, a =a , b =b, varx =var(reputationalcontrol.subset$reputationalframe), vary=var(reputationalcontrol.subset$petition1), 
             varm =var(reputationalcontrol.subset$policyapproval1), alpha = .05, interval = NULL)
####################Figure 15a############
ilcontrol.subset<-subset(data,ilframe==1|noframe==1)
model.m.ilcontrol<- lm(policyapproval1~ilframe, data = ilcontrol.subset)
model.y.ilcontrol<-glm(protest1~ilframe+policyapproval1, data = ilcontrol.subset, family=binomial(link="probit"))
mediatedilcontrol.protest<- mediate(model.m.ilcontrol, model.y.ilcontrol, treat = "ilframe", mediator = "policyapproval1", boot=TRUE, sims=2000,conf.level = .95)
summary(mediatedilcontrol.protest);plot(mediatedilcontrol.protest)

#il group vs. moral
ilmoral.subset<-subset(data,ilframe==1|moralframe==1)
model.m.ilmoral<- lm(policyapproval1~ilframe, data = ilmoral.subset)
model.y.ilmoral<-glm(protest1~ilframe+policyapproval1,  data = ilmoral.subset, family=binomial(link="probit"))
mediatedilmoral.protest<- mediate(model.m.ilmoral, model.y.ilmoral, treat = "ilframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedilmoral.protest); plot(mediatedilmoral.protest)

#il group vs. reputational
ilreputational.subset<-subset(data,ilframe==1|reputationalframe==1)
model.m.ilreputational<- lm(policyapproval1~ilframe, data = ilreputational.subset)
model.y.ilreputational<-glm(protest1~ilframe+policyapproval1, data = ilreputational.subset,family=binomial(link="probit"))
mediatedilreputational.protest<- mediate(model.m.ilreputational, model.y.ilreputational, treat = "ilframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedilreputational.protest); plot(mediatedilreputational.protest)

#moral group vs. control
moralcontrol.subset<-subset(data,moralframe==1|noframe==1)
model.m.moralcontrol<- lm(policyapproval1~moralframe, data = moralcontrol.subset)
model.y.moralcontrol<-glm(protest1~moralframe+policyapproval1,  data = moralcontrol.subset,family=binomial(link="probit"))
mediatedmoralcontrol.protest<- mediate(model.m.moralcontrol, model.y.moralcontrol, treat = "moralframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedmoralcontrol.protest);plot(mediatedmoralcontrol.protest)

#moral group vs. il group: no estimation -- this is just the inverse of il group vs. moral (above)

#moral group vs. reputational
moralreputational.subset<-subset(data,moralframe==1|reputationalframe==1)
model.m.moralreputational<- lm(policyapproval1~moralframe, data = moralreputational.subset)
model.y.moralreputational<-glm(protest1~moralframe+policyapproval1,  data = moralreputational.subset,family=binomial(link="probit"))
mediatedmoralreputational.protest<- mediate(model.m.moralreputational, model.y.moralreputational, treat = "moralframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedmoralreputational.protest); summary(mediatedmoralreputational.protest)

#reputational group vs. control
reputationalcontrol.subset<-subset(data,reputationalframe==1|noframe==1)
model.m.reputationalcontrol<- lm(policyapproval1~reputationalframe, data = reputationalcontrol.subset)
model.y.reputationalcontrol<-glm(protest1~reputationalframe+policyapproval1,  data = reputationalcontrol.subset,family=binomial(link="probit"))
mediatedreputationalcontrol.protest<- mediate(model.m.reputationalcontrol, model.y.reputationalcontrol, treat = "reputationalframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedreputationalcontrol.protest);plot(mediatedreputationalcontrol.petition)

#reputational group vs il group: no estimation -- this is just the inverse of il group vs. reputational (above)
#reputational group vs. moral group: no estimation -- this is just the inverse of moral group vs. reputational (above)
###Table 16a
sigma.x.ilcontrol<-sd(ilcontrol.subset$ilframe)#calculate std dev of x for power analysis
sigma.m.ilcontrol<-sd(ilcontrol.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.ilcontrol<-sd(ilcontrol.subset$protest1)#calculate std dev of y for power analysis
summary(n.il+n.control)#sample size to report
a<- summary(model.m.ilcontrol)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.ilcontrol)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.il+n.control), power = NULL, a =a , b =b, varx =var(ilcontrol.subset$ilframe), vary=var(ilcontrol.subset$protest1), 
             varm =var(ilcontrol.subset$policyapproval1), alpha = .05, interval = NULL)

#il-moral
sigma.x.ilmoral<-sd(ilmoral.subset$ilframe)#calculate std dev of x for power analysis
sigma.m.ilmoral<-sd(ilmoral.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.ilmoral<-sd(ilmoral.subset$protest1)#calculate std dev of y for power analyais
summary(n.il+n.moral)#sample size to report
a<- summary(model.m.ilmoral)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.ilmoral)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.il+n.moral), power = NULL, a =a , b =b, varx =var(ilmoral.subset$ilframe), vary=var(ilmoral.subset$protest1), 
             varm =var(ilmoral.subset$policyapproval1), alpha = .05, interval = NULL)

#il-reputational
sigma.x.ilreputational<-sd(ilreputational.subset$ilframe)#calculate std dev of x for power analysis
sigma.m.ilreputational<-sd(ilreputational.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.ilreputational<-sd(ilreputational.subset$protest1)#calculate std dev of y for power analyais
summary(n.il+n.reputational)#sample size to report
a<- summary(model.m.ilreputational)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.ilreputational)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.il+n.reputational), power = NULL, a =a , b =b, varx =var(ilreputational.subset$ilframe), vary=var(ilreputational.subset$protest1), 
             varm =var(ilreputational.subset$policyapproval1), alpha = .05, interval = NULL)

#moral-control
sigma.x.moralcontrol<-sd(moralcontrol.subset$moralframe)#calculate std dev of x for power analysis
sigma.m.moralcontrol<-sd(moralcontrol.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.moralcontrol<-sd(moralcontrol.subset$protest1)#calculate std dev of y for power analyais
summary(n.moral+n.control)#sample size to report
a<- summary(model.m.moralcontrol)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.moralcontrol)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.moral+n.control), power = NULL, a =a , b =b, varx =var(moralcontrol.subset$moralframe), vary=var(moralcontrol.subset$protest1), 
             varm =var(moralcontrol.subset$policyapproval1), alpha = .05, interval = NULL)

###moral-reputational
sigma.x.moralreputational<-sd(moralreputational.subset$moralframe)#calculate std dev of x for power analysis
sigma.m.moralreputational<-sd(moralreputational.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.moralreputational<-sd(moralreputational.subset$protest1)#calculate std dev of y for power analyais
summary(n.moral+n.reputational)#sample size to report
a<- summary(model.m.moralreputational)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.moralreputational)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.moral+n.reputational), power = NULL, a =a , b =b, varx =var(moralreputational.subset$moralframe), vary=var(moralreputational.subset$protest1), 
             varm =var(moralreputational.subset$policyapproval1), alpha = .05, interval = NULL)

##reputational-control
sigma.x.reputationalcontrol<-sd(reputationalcontrol.subset$reputationalframe)#calculate std dev of x for power analysis
sigma.m.reputationalcontrol<-sd(reputationalcontrol.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.reputationalcontrol<-sd(reputationalcontrol.subset$protest1)#calculate std dev of y for power analyais
summary(n.reputational+n.control)#sample size to report
a<- summary(model.m.reputationalcontrol)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.reputationalcontrol)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.reputational+n.control), power = NULL, a =a , b =b, varx =var(reputationalcontrol.subset$reputationalframe), vary=var(reputationalcontrol.subset$protest1), 
             varm =var(reputationalcontrol.subset$policyapproval1), alpha = .05, interval = NULL)

####################Table 17a############
ilcontrol.subset<-subset(data,ilframe==1|noframe==1)
model.m.ilcontrol<- lm(policyapproval1~ilframe, data = ilcontrol.subset)
model.y.ilcontrol<-glm(donate1~ilframe+policyapproval1, data = ilcontrol.subset, family=binomial(link="probit"))
mediatedilcontrol.donate<- mediate(model.m.ilcontrol, model.y.ilcontrol, treat = "ilframe", mediator = "policyapproval1", boot=TRUE, sims=2000,conf.level = .95)
summary(mediatedilcontrol.donate);plot(mediatedilcontrol.donate)

#il group vs. moral
ilmoral.subset<-subset(data,ilframe==1|moralframe==1)
model.m.ilmoral<- lm(policyapproval1~ilframe, data = ilmoral.subset)
model.y.ilmoral<-glm(donate1~ilframe+policyapproval1,  data = ilmoral.subset, family=binomial(link="probit"))
mediatedilmoral.donate<- mediate(model.m.ilmoral, model.y.ilmoral, treat = "ilframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedilmoral.donate); plot(mediatedilmoral.donate)

#il group vs. reputational
ilreputational.subset<-subset(data,ilframe==1|reputationalframe==1)
model.m.ilreputational<- lm(policyapproval1~ilframe, data = ilreputational.subset)
model.y.ilreputational<-glm(donate1~ilframe+policyapproval1, data = ilreputational.subset,family=binomial(link="probit"))
mediatedilreputational.donate<- mediate(model.m.ilreputational, model.y.ilreputational, treat = "ilframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedilreputational.donate); plot(mediatedilreputational.donate)

#moral group vs. control
moralcontrol.subset<-subset(data,moralframe==1|noframe==1)
model.m.moralcontrol<- lm(policyapproval1~moralframe, data = moralcontrol.subset)
model.y.moralcontrol<-glm(donate1~moralframe+policyapproval1,  data = moralcontrol.subset,family=binomial(link="probit"))
mediatedmoralcontrol.donate<- mediate(model.m.moralcontrol, model.y.moralcontrol, treat = "moralframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedmoralcontrol.donate);plot(mediatedmoralcontrol.donate)

#moral group vs. il group: no estimation -- this is just the inverse of il group vs. moral (above)

#moral group vs. reputational
moralreputational.subset<-subset(data,moralframe==1|reputationalframe==1)
model.m.moralreputational<- lm(policyapproval1~moralframe, data = moralreputational.subset)
model.y.moralreputational<-glm(donate1~moralframe+policyapproval1,  data = moralreputational.subset,family=binomial(link="probit"))
mediatedmoralreputational.donate<- mediate(model.m.moralreputational, model.y.moralreputational, treat = "moralframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedmoralreputational.donate); summary(mediatedmoralreputational.donate)

#reputational group vs. control
reputationalcontrol.subset<-subset(data,reputationalframe==1|noframe==1)
model.m.reputationalcontrol<- lm(policyapproval1~reputationalframe, data = reputationalcontrol.subset)
model.y.reputationalcontrol<-glm(donate1~reputationalframe+policyapproval1,  data = reputationalcontrol.subset,family=binomial(link="probit"))
mediatedreputationalcontrol.donate<- mediate(model.m.reputationalcontrol, model.y.reputationalcontrol, treat = "reputationalframe", mediator = "policyapproval1", boot=TRUE, sims=1000)
summary(mediatedreputationalcontrol.donate);plot(mediatedreputationalcontrol.petition)

#reputational group vs il group: no estimation -- this is just the inverse of il group vs. reputational (above)
#reputational group vs. moral group: no estimation -- this is just the inverse of moral group vs. reputational (above)
##########Table 18a########
sigma.x.ilcontrol<-sd(ilcontrol.subset$ilframe)#calculate std dev of x for power analysis
sigma.m.ilcontrol<-sd(ilcontrol.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.ilcontrol<-sd(ilcontrol.subset$donate1)#calculate std dev of y for power analysis
summary(n.il+n.control)#sample size to report
a<- summary(model.m.ilcontrol)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.ilcontrol)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.il+n.control), power = NULL, a =a , b =b, varx =var(ilcontrol.subset$ilframe), vary=var(ilcontrol.subset$donate1), 
             varm =var(ilcontrol.subset$policyapproval1), alpha = .05, interval = NULL)

#il-moral
sigma.x.ilmoral<-sd(ilmoral.subset$ilframe)#calculate std dev of x for power analysis
sigma.m.ilmoral<-sd(ilmoral.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.ilmoral<-sd(ilmoral.subset$donate1)#calculate std dev of y for power analyais
summary(n.il+n.moral)#sample size to report
a<- summary(model.m.ilmoral)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.ilmoral)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.il+n.moral), power = NULL, a =a , b =b, varx =var(ilmoral.subset$ilframe), vary=var(ilmoral.subset$donate1), 
             varm =var(ilmoral.subset$policyapproval1), alpha = .05, interval = NULL)

#il-reputational
sigma.x.ilreputational<-sd(ilreputational.subset$ilframe)#calculate std dev of x for power analysis
sigma.m.ilreputational<-sd(ilreputational.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.ilreputational<-sd(ilreputational.subset$donate1)#calculate std dev of y for power analyais
summary(n.il+n.reputational)#sample size to report
a<- summary(model.m.ilreputational)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.ilreputational)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.il+n.reputational), power = NULL, a =a , b =b, varx =var(ilreputational.subset$ilframe), vary=var(ilreputational.subset$donate1), 
             varm =var(ilreputational.subset$policyapproval1), alpha = .05, interval = NULL)

#moral-control
sigma.x.moralcontrol<-sd(moralcontrol.subset$moralframe)#calculate std dev of x for power analysis
sigma.m.moralcontrol<-sd(moralcontrol.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.moralcontrol<-sd(moralcontrol.subset$donate1)#calculate std dev of y for power analyais
summary(n.moral+n.control)#sample size to report
a<- summary(model.m.moralcontrol)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.moralcontrol)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.moral+n.control), power = NULL, a =a , b =b, varx =var(moralcontrol.subset$moralframe), vary=var(moralcontrol.subset$donate1), 
             varm =var(moralcontrol.subset$policyapproval1), alpha = .05, interval = NULL)

###moral-reputational
sigma.x.moralreputational<-sd(moralreputational.subset$moralframe)#calculate std dev of x for power analysis
sigma.m.moralreputational<-sd(moralreputational.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.moralreputational<-sd(moralreputational.subset$donate1)#calculate std dev of y for power analyais
summary(n.moral+n.reputational)#sample size to report
a<- summary(model.m.moralreputational)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.moralreputational)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.moral+n.reputational), power = NULL, a =a , b =b, varx =var(moralreputational.subset$moralframe), vary=var(moralreputational.subset$donate1), 
             varm =var(moralreputational.subset$policyapproval1), alpha = .05, interval = NULL)

##reputational-control
sigma.x.reputationalcontrol<-sd(reputationalcontrol.subset$reputationalframe)#calculate std dev of x for power analysis
sigma.m.reputationalcontrol<-sd(reputationalcontrol.subset$policyapproval1)#calculate std dev of mediator for power analysis
sigma.y.reputationalcontrol<-sd(reputationalcontrol.subset$donate1)#calculate std dev of y for power analyais
summary(n.reputational+n.control)#sample size to report
a<- summary(model.m.reputationalcontrol)$coefficients[2] # path between predictor and mediator
b<- summary(model.y.reputationalcontrol)$coefficients[3] #path between mediator and dependent variable
wp.mediation(n = (n.reputational+n.control), power = NULL, a =a , b =b, varx =var(reputationalcontrol.subset$reputationalframe), vary=var(reputationalcontrol.subset$donate1), 
             varm =var(reputationalcontrol.subset$policyapproval1), alpha = .05, interval = NULL)



